arXiv: 1504.03158vl [quant-ph] 13 Apr 2015 


Sued et al. 


RESEARCH 


Quantum Lattice Boltzmann is a quantum walk 

Sauro Sued 1 *, Francois Fillion-Gourdeau and Silvia Palpacelli 2 


Correspondence: 
succi@iac.rm.cnr.it 
1 IAC-CNR, Istituto per le 
Applicazioni del Calcolo “Mauro 
Picone” , Via dei Taurini, 19, 
00185 Rome, Italy 
Full list of author information is 
available at the end of the article 


Abstract 

Numerical methods for the 1-D Dirac equation based on operator splitting and 
on the quantum lattice Boltzmann (QLB) schemes are reviewed. It is shown that 
these discretizations fall within the class of quantum walks, i.e. discrete maps for 
complex fields, whose continuum limit delivers Dirac-like relativistic quantum 
wave equations. The correspondence between the quantum walk dynamics and 
these numerical schemes is given explicitly, allowing a connection between 
quantum computations, numerical analysis and lattice Boltzmann methods. The 
QLB method is then extended to the Dirac equation in curved spaces and it is 
demonstrated that the quantum walk structure is preserved. Finally, it is argued 
that the existence of this link between the discretized Dirac equation and 
quantum walks may be employed to simulate relativistic quantum dynamics on 
quantum computers. 
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1 Introduction 

A quantum walk (QW) is defined as the quantum analogue of a classical random 
walk, where the “quantum walker” is in a superposition of states instead of being 
described by a probability distribution. One of the earliest realization of this concept 
was proposed by Feynman as a discrete version of the massive Dirac equation [1]. In 
recent times, there has been a surge of interest for this topic, due to the conceptual 
and possible practical import of QW’s as discrete realizations of stochastic quantum 
processes and because they can solve certain problems with an exponential speedup, 
i.e. using exponentially less operations than classical computations [2]. Moreover, 
QW’s are amenable to a number of experimental realizations, such as ion traps 
[3, 4, 5], liquid-state nuclear-magnetic-resonance quantum-information processor 
[6], photonic devices [7] and other types of optical devices [8]. As a result, they 
hold promise of playing an important role in many areas of modern physics and 
quantum technology, such as quantum computing, foundational quantum mechanics 
and biophysics [9]. 

One of the most interesting features of discrete QW’s is their continuum limit, 
which recovers a broad variety of relativistic quantum wave equations [10, 11, 12]. As 
stated earlier, this was first discussed by Feynman and is now known as Feynman’s 
checkerboard [13]. This was originally formulated for the free Dirac equation but 
extensions of these ideas, which include the coupling to external fields, have been 
investigated [10, 11, 12]. Pioneering work have been performed in which the Dirac 
equation is related to cellular automata [14, 15]. Lately, the link between QW’s and 
the Dirac equation have been discussed extensively [10, 11, 12, 16]. In these studies, 
the starting point is a general QW formulation from which the continuum limit is 
evaluated and then related to the Dirac equation. 
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2 QUANTUM WALKS 


In this article, QW’s and their relations to some known numerical schemes of the 
Dirac equation are reviewed from a slightly different perspective: it will be demon¬ 
strated that the most general QW’s are obtained from lattice discretizations of 
the relativistic quantum wave equation for spin-1/2 particles. More precisely, start¬ 
ing from the continuum Dirac equation, it is shown that QW’s can be placed in 
one-to-one correspondence with numerical schemes based on operator splitting and 
the QLB scheme. These numerical methods have been developed and employed as 
efficient numerical tools to solve relativistic quantum mechanics problems on clas¬ 
sical computers [17, 18, 19]. They have a number of interesting properties: they are 
easy to code, they can be easily parallellized and are very versatile. Moreover, their 
mathematical structure and the fact that the time discretization is realized by a 
set of unitary transformations makes the link to QW’s possible. This connection is 
explored below and represents one of the main purposes of this article. Much of the 
material discussed here, in particular the numerical simulations, is not completely 
new, which is in line with the main purpose of this paper, namely an attempt to 
bridge the techniques utilized in numerical analysis and (quantum) Lattice Boltz¬ 
mann theory to the language of quantum computing. 

This article is organized as follows. In Section 2, a general formulation of QW 
is presented, where the transfer matrix is time and space dependent. In Section 3, 
the split operator method for the Dirac equation is presented, along with its exact 
correspondence with QW. Section 4 is devoted to the QLB method and connections 
with QW. Section 5 is devoted to a qualitative discussion of the link between these 
numerical schemes and quantum computation. In Section 6, the schemes are casted 
in the form of a propagation-relaxation process and the notion of quantum equi¬ 
librium is introduced. Based on the analogy between QW and QLB, a new QLB 
scheme for the (1+1) Dirac equation in curved space is proposed in Section 7. Fi¬ 
nally, the generalization of these methods to many dimensions is briefly discussed 
and numerical results are presented in Section 8 

2 Quantum walks 

Let us consider a (1 + 1) quantum walk on the line for a pair of complex wave 
functions (bi-spinor ?/>), obeying a discrete space-time evolution equations described 
by the following discrete map [10, 11, 12, 20]: 


kr] 

= B 

+-1 

krj 


7+1. 


Here, the indices j, ra G N 0 Z label points on a discretization of space and time, 
respectively. The object Bj, n is a two-by-two matrix with components [12] 


e l0i i’ n cos 6j^ n e lf3 ^ n sin <9y n 
—e~ % Pj' n sin 0j :Tl e~ ta *> n cos #y n 


( 2 ) 


This matrix is a U( 2) operator (B G 17(2)) parametrized by the three space-time 
dependent Euler angles 0j, n ,aj, n ,(3j, n and a space-time dependent phase £y n . The 
latter is relevant when it depends on time and space, i.e. when it is local. If it is 
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global (Cj,n = £), it disappears from any observables and becomes unimportant. 
This occurs because when the phase is space-time dependent, it does not commute 
with the time and space translation operators. 

The matrix obeys B E SU (2) only when = kir for all j, n, with k E Z. Thus, 
this formulation is slightly more general than QW’s considered in [10, 11, 20] where 
B E SU( 2) is studied. As seen in the next section, the choice B E U(2) will be 
important to have a general connection between mass terms of the Dirac equation 
and QW’s. Finally, the U( 2) QW can also be implemented on quantum computers 
because the matrix B is a unitary transformation: it represents the most general 
QW consistent with quantum computations. 

In the above, the amplitudes Vh ,2 code for the probability of the quantum walker 
to move up (down) along the lattice site j E Z at the time step n E N. This 
is a very rich structure, which has been shown to recover a variety of important 
quantum wave equations, as soon as the Euler angles are allowed to acquire a space- 
time dependence [12]. In addition, it provides a wealth of potential algorithms for 
quantum computing. This was studied extensively in [10, 11, 12] by analysing the 
continuum limit of these QW, yielding different versions of the Dirac equation. In 
this work, the opposite path is taken: it is shown that specific discretizations of 
the Dirac equation, using either a split-operator approach or the lattice Boltzmann 
technique, naturally lead to a QW formulation. 

3 Split-operator and quantum walks 

The starting point of this discussion is the 1-D Dirac equation in Majorana repre¬ 
sentation written as (in units where c = h = 1): 

idt'ip^z, t) = [ ~icr z d z + M(z, t)] t/>(z, *), (3) 

with the bi-spinor ^ E Z/ 2 (M, C 2 ). The generalized “mass” matrix M is space and 
time dependent and may include contributions from the physical mass, the coupling 
to an electromagnetic potential or any other type of coupling. One requirement, 
however, is that M is a Hermitian local operator without any derivatives. Generally, 
it can be written as 


M(z,t) = hM 0 (z, t) + a • M, 

= hM 0 (z,t) + a x M x (z,t) + a y M y (z,t ) + a z M z (z,t), 


( 4 ) 

( 5 ) 


where I 2 is the two-by-two identity matrix, (< Ji)i= x , y , z are Pauli matrices and the 
coefficients Mo )X) y )2 represent the time- and space-dependent external fields which 
couple to the spinor. 

The formal solution of Eq. (3) is given by 


t) = T exp 


—A ta z d z —i f M(z,t')dt' 
Jt 0 




(6) 


where T is the time-ordering operator, to is the initial time and At = t — to. Using 
an operator splitting technique, the solution in Eq. (6) can be approximated by [19] 


/ i/j( z , t) = exp [— iAtM(z, to)] exp [— Ata z d z ] i/j(z, to) + 0(At 2 ). 


( 7 ) 
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3 SPLIT-OPERATOR AND QUANTUM WALKS 


The first exponential is a translation (streaming) operator which shifts the spinor 
components according to: 


exp [—Atcr z d z ] ip(z, to) 


fpi 0 - A t,t 0 ) 
ip 2 (z + At,t 0 ) 


(8) 


This suggests to use a spatial discretization where Az = At (this corresponds to a 
Courant-Friedrichs-Lewy (CFL) condition C = cAt/Az = 1, c being the speed of 
light) such that the translation is exact on the lattice. Eq. (7) can then be written 
as: 


,n+1 
1 ,3 
,n+1 
2,j . 


exp [— iAtMj^ n \ 


Mj-i 


( 9 ) 


where we defined the following quantities on the lattice: M^ n := M(nAt,jAz) 
and x/j'j := ^(nAt, jAz). This last equation yields a numerical scheme to solve 
the Dirac equation. This numerical scheme has interesting properties, as discussed 
extensively in [19, 17, 21]: it can be extended to higher order accuracy, it can be 
easily parallellized and it can be easily coded on a computer. In the following, it is 
also demonstrated that it is completely equivalent to the U (2) QW described in the 
last section. 

Eq. (9) is in the form of Eq. (1). Moreover, the exponential is also a unitary 
matrix: 


B' := exp [— iAtMj, n ] G U( 2), (10) 

and thus, there clearly exists a connection between B and B '. They are expressed 
in different representation: B uses the Euler angle parametrization while B' is ex¬ 
pressed in the canonical representation obtained by the exponential mapping of the 
Lie algebra. The latter is given explicitly by 


B' 


exp [—iAtMoj jn ] exp [—iAta • My n ] 


exp [~iAtM 0 j, n ] 


h cos(|Mj ?n |At) 


(T • M i>n 

|-M-2,ra | 


sin(|Mj >n |At) 


( 11 ) 

( 12 ) 


where 


|M 


■3,n 


= ^ M Ln 


■M 2 - 

V,3,n 


M 2 ■ . 

z,j,n 


(13) 


It is well-known that parametrizations of U( 2) matrices are related to each other 
[22] and it can be determined that the mapping between both representations is 
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given by 


£j,n 

II 

i 

3 

> 

(14) 

tan (a j>n ) 


(15) 

tan (fan) 

~ M ■ ’ 

lvl y,j,n 

(16) 

tan (9j >n ) 


(17) 

When Mq = 0, 

one recovers the SU{2 ) QW because then, B := B'\m 0 =o 

G S'(7(2). 


These last equations give a one-to-one correspondence between the QW formula¬ 
tion, characterized by the parameters oy, n , and the discretization of 

the Dirac equation with a generalized mass term Mj ?n . Therefore, the last results 
show that in the continuum limit, every space-time dependent QW on the line be¬ 
comes a time-dependent ID Dirac equation with a specific mass matrix. A given QW 
can thus be fully characterized by the relativistic dynamics of an electron coupled 
to a space-time dependent external field. This occurs because there is an equiva¬ 
lence between the discretization of the Dirac equation, based on operator splitting, 
and the QW formulation. Moreover, this connection may be the base for the im¬ 
plementation of a quantum algorithm that solves the Dirac equation on quantum 
computers. 

The operator splitting technique presented here also bears a close relationship 
with the QLB technique, to which we now turn. 

4 Quantum lattice Boltzmann, operator splitting and Quantum 
Walks 

The QLB was inspired by a direct analogy between the way the Dirac equation goes 
to the Schroedinger equation in the limit v/c —> 0, and the way that the Navier- 
Stokes equations of classical fluid-dynamics emerge from the Boltzmann equation 
in the limit of small Knudsen number, Kn —>• 0, where Kn = l/L is the ratio 
of the molecular mean free path to the typical macroscopic scale. In both cases, 
the smallness parameter controls the enslaving of the fast modes to the slow ones: 
non-equilibrium to equilibrium for the classical case, versus excited states to ground 
state in the quantum one. Of course, the quantum case shows no genuine relaxation 
since its dynamics is reversible. Yet, enslaving can be interpreted in the sense of fast 
oscillations around a local quantum equilibrium (Zitterbewegung), which average 
out once time is coarse-grained on a scale larger than the period of the fast oscilla¬ 
tions. So, Zitterbewegung may be regarded as the quantum relativistic analogue of 
classical non-equilibrium fluctuations. 

Based on this analogy, QLB was formulated as a lattice Boltzmann analogue of 
the Dirac equation in the Majorana representation, where the streaming matrix is 
real. To obtain the QLB scheme, it is convenient to write the Dirac equation as 

[d t +V a d z ] ipa(z,t) = -iJ2M ab (z,t)Mz,t), 

b 


(18) 
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where a, b are spinor indices and the “microscopic velocities” are given by v\^ = 
d=l. This is clearly in a “Boltzmann-like” form with two discrete velocities and a 
collision term M. When QLB is employed in fluid mechanics to solve the continuum 
equations of motion, there is a family of possible numbers of discrete velocities for a 
given lattice [23] and each choice yields a different numerical scheme (for instance, 
the 9 velocities scheme in 2-D and the 27 velocities scheme in 3-D are popular 
choices on square lattices [24]). For the Dirac equation, this choice is dictated by 
the mathematical structure of the equation. 

The lattice Boltzmann equation is then written as [25, 26] 

^i{z + v a At, t + At) = (z, t) - iAt ^2 M ab(z , t)^ b (z, t ) + 0(At 2 ). (19) 

b 

As for the splitting method described in the last section, this suggests to use a space 
discretization where Az = At. Using a “naive” approach in line with LB methods 
in fluid mechanics, the last equation would become 


U n+1 

™i4+i 

? /; n + 1 
/ 2,^ — 1 


A"; 


— iAtMj^ n 




( 20 ) 


Then, the matrix —iM on the right acts as collision operator while the streaming 
is executed by the left part of the equation. This scheme is derived by using the 
formal analogy between the Dirac equation, the Boltzmann equation and the LB 
technique. However, the resulting numerical method is unstable and the L 2 norm 
is not preserved [27]. 

It is possible, however, to recover a stable and norm-preserving method. Instead 
of using the “naive” LB equation discretization given in Eq. (20), the mass term on 
the right-hand side of Eq. (19) is discretized by using an implicit Crank-Nicolson 
average: 


iM(z, t) 


1pl(z,t) 

(Z,t) 



ib n+1 

™14+1 

q /; n + 1 

A 2 4-i 


r hj ] \ 

Mj\ f' 


( 21 ) 


Reporting this into Eq. (19), one obtains the second order accurate QLB scheme: 


\b n+1 

— T 

"< 4 / 

1 

— J-j,n 



+ 0(At 2 ), 


( 22 ) 


where the transfer matrix is given by 



At 

-1 

At 

T — 

J-j.n — 

I 2 + i-prMj 

2 


I 2 

2 


To link these results with the ones of the last sections, it is convenient to shift the 
spinor components by one lattice point (we let ^ 14+1 —^ Vh 4 and ^ 24-1 —^ ^ 24 ). 
Then, the QLB scheme becomes 


<44 

.<44 


= T 


3,n 


<4i-i 


+ 0(At 2 ), 


(24) 


























4.1 An explicit example: the free case 
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which is in the same form as the QW and the splitting method presented previously. 
The transfer matrix can be evaluated explicitly. It is given by 


Tj ’ n ~ Ci 


1 - iAtM zJ>n + ^M| n 


At [iM x ,j,n T MyJ i n') 
1 + iAtM zJt n + 


(25) 


where C hJl = 1 + iAtM( ) j n — M 2 n and M 2 n = Mqj u — Mj n . It can be readily 
verified that T^ n Tj n = I 2 . As a consequence, the transfer matrix is unitary T ] n G 
17(2), for any size of the time step and thus, conserves the L 2 norm. Moreover, just 
like the splitting method, there exists a correspondence with QW’s because both 
have U(2) collision matrices. The identification with the general QW in Eq. (1) 
yields 


tan(£j,n) 


tan(afy n ) 
tan (/3 i>n ) 


tan (0j >n ) 


Im (Cj >n ) AtM 0 j, r . 


R e(C j>n ) 1 - fM^ n 


AtM, 


z,3,n 


1 + A£m 2 ’ 

1 4 J,n 


Mx 


M 


y,3,n 


-At 


M 2 • + M 2 • 


N ( 1 + ^f M ln ) 2 + At2M Ln 


(26) 

(27) 

(28) 
(29) 


These relations give a correspondence between the QLB technique and the QW. 
They are similar to the ones for the splitting method displayed in Eqs. (14) to (17) 
and actually serve the same purpose: they allow to map a numerical scheme to 
the QW. The differences are obviously due to the fact that QLB is based on a LB 
formulation combined with a Crank-Nicolson average to insure stability while the 
splitting scheme separates the Dirac equation into different operators which can be 
integrated exactly. However, both methods share the same general structure where 
a streaming step is followed by a collision step. 


4.1 An explicit example: the free case 

The Dirac equation for the massive and free case (no external field coupling) is 

id t ip(z, t ) = [-id z d z + (j y m\ ^(z, t), (30) 

where M y = m is a physical fermion mass. This representation of the Dirac equation 
yields real spinor components, as readily seen by writing the equation component¬ 
wise: 


(d t + d z )ipi(z, t) = —ra^OM), 
(d t -d z )^ 2 (z,t) = +rra/q(£,£). 


(31) 

(32) 
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5 PROSPECTS FOR QUANTUM SIMULATION 


The QLB scheme, for the massive and free case, reads as follows [18]: 




= T{„ 


^■44 


(33) 


where the transfer matrix is can be obtained from Eq. (25) by setting Mq = M x = 
M z = 0 and M y = m: 


Tf re e 


a(m) b(m) 
—b(m) a(m) 


1 

-| | A t 2 m 2 

1 ' 4 


A t 2 m 2 


mAt 


—rriAt 
A t 2 m 2 


(34) 


It is readily seen that a(m) and b(m) are second-order Pade’-like approximants of 
cos(ra) and sin(ra), respectively. This is the natural consequence of the implicit 
time-marching (Crank-Nicolson) scheme as applied to the (1+1) Dirac equation in 
Majorana form. 

For a free particle, the correspondence to QW is particularly simple and is given 
by 


(35) 

(36) 


This mapping of the QLB to the QW is exact for any value of m. 


£ = a = f3 = 0 , 

_ mAt 

tan(0) = 


1 - 


A t 2 r 


5 Prospects for quantum simulation 

The numerical methods described in Sections 3-4 can be straightforwardly imple¬ 
mented on classical computers. However, the stream-collide structure of these nu¬ 
merical scheme makes them suitable for an efficient implementation on quantum 
computers as well. In particular, they can be written as: 

y n+1 = -Vyr''. (37) 


where Sj is the shift operator, defined as: 


Srfj 


+i-l 
+i+1. 


(38) 


This shift operator is a unitary operation and it can be realized experimentally 
by using fundamental quantum gates [6, 20]. The rotation operator Bj jU belongs 
to 1/(2) and therefore can also be realized by these quantum gates, as any other 
unitary transformations [28]. Therefore, it is possible to map the numerical method 
to quantum walk, which can be implemented efficiently on quantum computers. This 
mapping is possible because each step of the scheme is a unitary transformation: this 
makes these schemes norm-preserving and sets the link with quantum computations. 
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The latter would be particularly useful for the study of relativistic quantum sys¬ 
tems where a time-dependent solution of the Dirac equation is required, such as in 
very high intensity laser physics [29] or graphene physics [30]. 

Another subject of major interest for future research is the extension of the 
QLB methodology to quantum many body systems and quantum field theory, two 
paramount sectors of modern physics which are particularly exposed to the limi¬ 
tations of classical (non-quantum) electronic computing. Progress in this direction 
depends on the ability to replace the quantum wavefunction by the corresponding 
second-quantized quantum operators, and show that the dynamics of the second- 
quantized QLB scheme still preserves the appropriate equal-time commutation re¬ 
lations. Preliminary efforts along this line have been developed in [31] in 1 + 1 
dimensions. Extensions to strongly non-linear field theories in d > 1 remain to 
be explored. As to quantum-many body problems, LB-like methods have been re¬ 
cently adapted to electronic structure simulations [32] In this work, a classical LB 
scheme is employed to solve the Kohn-Sham equations of density functional theory 
in the form diffusion-reaction equations in imaginary time. Allied QLB schemes 
could prove very useful to solve the corresponding real-time quantum many-body 
transport problems within the framework of time-dependent density functional the¬ 
ory. 

Finally, we wish to point out the intriguing possibility of realizing both quantum 
and classical LB schemes on quantum analogue simulators, as recently explored in 

[33]. 

6 Quantum equilibria 

In view of quantum computing implementations, it is of interest to cast the Dirac 
equation in the form of a propagation-relaxation process, where the collision matrix 
is now interpreted as a scattering process, relaxing the spinor component around a 
local quantum equilibrium. 

For this purpose, it is useful to reconsider the 1-D Dirac equation in the form: 

[d t + v a d z } ip a (z , t) = -i^2 M ab (z, t)ip b (z, t). (39) 

b 

Then, it is formally possible to define a local equilibrium as: 

ip eq (z,t) :=Uip(z,t), (40) 

where U is a unitary matrix that depends on M. This transformation is chosen to 
recast the Dirac equation in relaxation form: 

[dt + cr z d z \ t ) = -fi t ) - 'Ipeqiz, t)] , (41) 

where Q = iM[I + C/] -1 . The explicit value of I/, Q is not unique but a convenient 
choice is U = e lMr . 

In this vests, the Dirac equation looks formally like a linear Boltzmann equation 
for two-component models in the single relaxation time approximation [25] . There¬ 
fore, it can be interpreted as a propagation-relaxation process in imaginary time, 
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6 QUANTUM EQUILIBRIA 


whereby collisions, implemented by the scattering operator drive oscillations 
around the equilibrium distribution ^eq- By defining a post-collision wavefunction 
as: 


ip (z, t ) := (1 — A tfl) ip(z , t) + A tLhp eq (z, £), (42) 

the Dirac equation takes the most compact form 

ip(z + v a At, t + At) = ip (z, t). (43) 

This is particularly suitable to quantum computing implementations in the form of 
a classical stream-collide dynamics. The collision (relaxation) gate transforms the 
pre-collisional spinor ip into the post-collisional state ip ', and the streaming gate 
moves the post-collisional spinor to its destination location z±Az. Both operations 
are unitary and can be encoded in logical gates for quantum computing purposes 
[34, 35]. 

The expression (40) shows that the local equilibria is a linear function of the 
actual wavefunction ip, hence itself a function of space and time. 

The question is: will the actual wavefunctions ever reach this moving target, i.e, 
^ = V’eq? 

Based on its definition, this can only occur once ip eq lies in the null-space of the 
scattering matrix M, namely: 

M^eq = 0. (44) 

For the case of a free massive particle, it can be checked explicitly that the only 
solution is the trivial vacuum ip eq = 0. 

This means that the spinorial wavefunction is a superposition of (both slow and 
fast) zero-average oscillations around a local equilibrium which, consistently with 
the reversible nature of quantum mechanics, is actually never attained. 

Although this remains to be checked in detail, we conjecture that the same holds 
true for the case of a massive particle in an external potential, because in this case 
the Dirac equation is still linear. 

Based on (44), the condition for the local equilibrium to depart from a trivial 
vacuum is that the matrix M be singular, i.e. the local equilibrium is a zero-mode 
of the scattering matrix. 

Non-trivial quantum zero-modes, V’eq ^ 0, may indeed arise for non-linear quan¬ 
tum wave equations, such as the Gross-Pitaevski or the mean-field version of the 
Nambu-Jona-Lasinio model, to be dicussed shortly. A non-trivial local equilibrium 
would then signal a spontaneously broken symmetry, which is indeed the distinctive 
trait of the aforementioned non-linear quantum wave equations. 

Even though the notion of quantum equilibrium remains purely formal in nature, 
it is argued that it might nonetheless facilitate quantum computing implementations 
based on the compact expressions (42) and (43). This stands out a very interesting 
topic for future research. 
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7 QLB in curved space-time 

Quantum walks have been shown to map into Dirac-like equations in curved space 
as well by evaluating the continuum limit of certain QW’s [36, 12]. Here, in the 
same spirit as other sections, a QW structure is obtained by discretizing the Dirac 
equation in curved space time using a QLB-like approach. However, because the 
wave function propagates on a curved manifold, the structure of the resulting QW 
is different from Eq. (1) and should include a residency matrix that corrects the 
streaming step, which is strictly valid only in flat space. This is different from the 
result obtained in [36, 12]. 

The Dirac equation in a static (1+1) curved space writes as: 

7 a K( d n - + ^fi ,_1/2 ^(fl 1/2 e^)V’)] = -irmp, (45) 

where a = 0,1 and fi = t,z. In the above is the two-dimensional vierbein 
(zweibein) relating the components of the locally tangent Minkowski space basis 
(eo, ei), described by the metric 77 ^, to the global space basis (e t , e z ), described by 
the general metric g^ v . Also, g := det(g /JjU ) is the determinant of the metric tensor. 

The general form of the corresponding partial differential equation is 

dt^ip + a z A(z)d z vjj = Q(z, t)^, (46) 

with A(z) G M a function of space and Q(z, t ) the two-by-two gravitational collision 
matrix associated with Eq. (45). This becomes a hyperbolic system of equations 
where the advection speed A(z) should not be confused with the vector electrody¬ 
namic potential. Since the advection term is heterogeneous, a strict QLB structure, 
i.e. streaming along constant lightcones Az = +cAt, is no longer viable. Here, the 
situation is very similar to classical LB schemes on non-uniform grids (unsurpris¬ 
ingly, since in dimension D = 1, gravity is basically a stretching of the metric). 
For this purpose, several finite-volume LB (FVLB) schemes have been formulated, 
whose main outcome is that the streaming operator is no longer a diagonal matrix 
with eigenvalues ±1, as required by the formal QLB structure [37, 38]. 

In this respect, a finite-volume QLB scheme for the Dirac equation in curved space 
can be obtained by writing the Dirac equation as 

dt^ip + d z [cr z A(z )t/>] = Q'{z, (47) 

with Q'(z,t) := Q(z,t) + d z A(z)a z . Then, this equation is integrated over a control 
volume V (enclosed by a surface S ), extending from (j — |)A z to (j + |)A z. Fi¬ 
nally, applying a Crank-Nicolson time-marching and combining with upwind finite- 
differences, the discretized equation is given by 

WP-Mj + lfAfadS = y c (Qnj- 1^-1+OiyC) 

+ ^( ( 2l2J+lV , 2,j+l + ( 3l2,i' i / , 2,t 1 )’ (48) 

Wf - Ms -\j Mids = L (Q^j+Mj+i + Qnjl’fi 1 ) 


( 49 ) 
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7 QLB IN CURVED SPACE-TIME 


where c := Az/At is the uniform lattice light speed (or CFL condition). The bound¬ 
ary integrals are given by 

j> Aipi i2 dS = [Aip i, 2 ]j + i - [Aipi^j-i, (50) 

an thus, require an interpolation from the cell centers j, to the north and south 
boundaries at j db 1/2, respectively. Following a common practice in finite-volume 
formulations of hyperbolic problems, the flux terms are approximated by [39]: 


[Mi\ j+ i = Mij, ( 51 ) 

[Ail> i]j_i = (52) 

= AriV£j+i, (53) 

[Aip 2 ]j -1 = Ajiplj. (54) 


As a result: 


1 - 


Q'n,j \ 

2c J 


m 1 


/ n +1 

2c ^ 


1 - 


Q' 2 


22,j 

2 c 


o/j n + 1 

V2 ,j 


Q' 2 


21,j 


2 c 


n-\-l 

3 





Q\ 


in 


2 c 


$ 


2J+1> 



. ^21,j + l 
2 c 




n 

1,3 + 1' 


2c 


Q22J-1 

2c 



^,-1 


( 55 ) 





n 

2,3-1 


(56) 


It is readily seen that this reduces to a standard QLB in the limit of a uniform 
grid, when Aj/c = 1. In this case the streaming is diagonal with speed ±c and the 
spinors at (j, n + 1) are connected to the corresponding spinors at ( j d= l,n) by a 
local 2x2 matrix, which can be readily inverted to deliver a fully explicit map. 
However, when Aj/c ^ 1, the spinors at (j, n) also enter the map, so that local 
inversion delivers a slightly more elaborated structure, namely: 


■0” +1 = (R + TS)ipJ, 

which can be written more explicitly as 


(57) 


kf] 

= R 

V’i j 

+ T 

1 

5-0. 

1 

-1 

W\ 


Nr 




(58) 


In the above T is the local 2x2 transfer matrix including collisions, S is the 
streaming operator and R the local residency matrix, expressing the fraction of 
spinors which are left in the cell centered about z as the quantum system advances 
from t to t + At. This is depicted in Fig. 1. 

Clearly, the residency matrix vanishes in the case of a uniform mesh, i.e. no grav¬ 
ity. The mapping (58) represents the “gravitational” QLB. The detailed expressions 
of the streaming and residency matrices depend on the specific form of the metric 
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Figure 1 Sketch of the transfer and residency matrix elements. T\\ is the fraction of up-moving 
spinor jumping from j — 1 at time n — 1 to j at time n, while Rn is the fraction left in j — 1. T 22 
and R 22 bear the same meaning for the down moving spinor. 


tensor and associated vierbeins. Moreover, this analysis concentrates on the mathe¬ 
matical structure (streaming and collision steps) of the resulting scheme rather than 
on its numerical properties (convergence, stability, etc). These topics shall make the 
object of a future publication. 

8 Multi-dimensions 

The discretization presented in this work extends to the D + 1 dimensional case 
by applying the notion of operator splitting. This implies the inclusion of a new 
dynamic step which is entirely quantum: namely a ’’rotation”, designed so as to 
keep the spin aligned with the momentum along each of the three spatial directions. 
Schemes using this strategy can be found in [18, 19] and in [17] for higher order 
splittings. 

It might be that such rotation is not needed by formulating the Dirac equation as 
a random walk on other lattice with more natural topologies (the diamond) lattice. 
The QLB-QW equivalence in multi-dimensions will be discussed in a future publi¬ 
cation. However, to demonstrate the strength of the numerical schemes presented 
here and to show some possible applications for quantum computing, numerical 
results in 2-D are presented in the following. 

8.1 Numerical results 

As an example of possible applications of QLB scheme, we present two representa¬ 
tive simulations: Klein tunnelling in the presence of random impurities and Dirac 
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8 MULTI-DIMENSIONS 


equation with Nambu-Jona-Lasinio (NJL) interactions in 2 + 1 space-time dimen¬ 
sions. Details on the numerical methods used to obtain these results are given in 
[18, 19, 40]. Also, these results are not completely new as similar systems have been 
studied in [19, 41]. 

8.1.1 Graphene with random impurities 

In the first numerical test, the propagation of a Gaussian wave packet through a 
graphene sample with randomly distributed impurities is simulated [41]. In Ref. 
[41], simulations are performed for different values of the impurity concentration 
and the potential barrier, in order to provide an estimate of the effect of impurity 
concentration on the conductivity of the graphene sample. In Fig. 2, we report 
some representative snapshots of the first 1800 time steps of the simulation, at an 
impurity percentage^ 0.5% and V = 50 MeV. A lattice of size 2048 x 512 cells is 
used and the cell size is chosen to be Az = 0.96 nm, while the spreading of the initial 
Gaussian wave packet is a = 48 (in numerical units), leading to a Fermi frequency 
kp = 0.117 (80 MeV in physical units). In this simulation, a fully relativistic particle 
(m = 0) is considered. 

From Fig. 2, we can see that the wave packet is scattered by the impurities, giving 
rise to a plane front out of the initial Gaussian configuration. As a consequence of 
the randomness induced in the wave function by the disordered media, there is a 
momentum loss and therefore the motion of the wave packet is found to experience 
a corresponding slow down. It is also found that the wave packet takes more time 
to regroup as the impurity concentration and impurity potential are increased. 

8.1.2 Nambu-Jona-Lasinio interaction 

As a second example, we present a 2 + 1 space-time simulation of the Dirac equation 
with a Nambu-Jona-Lasinio (NJL) interaction [42]. The Dirac equation with a NJL 
interaction term driven by the coupling parameter g , reads as follows 

(i^d^ - m)ijj + + (V^VOtV) = 0, (59) 

where ijj = (Vh 5 ' 02 , V*3j ^a) T , 7 5 = ^7°7 1 7 2 7 3 and ijj = 7 /^ 7 °. 

This model represents a paradigm for dynamic mass acquisition via spontaneous 
symmetry breaking due to the non-linear interactions. 

Let us consider an initial condition given by the following Gaussian minimum- 
uncertainty wave packet: 

G 0 (z,y ) = (27r<7 2 )~ 1/2 exp(- Z ), (60) 

centered about (z,y) = (0,0), with initial width a. Let k z and k y be the initial 
energy of the wave packet and impose the following initial condition: 

ui(z,y) = u 2 (z,y) = C u G 0 (z,y) exp(i(k z z + k y y)), 

(61) 

di(z,y) = d 2 (z,y) = C d G 0 (z,y) exp(-i(k z z + k v y)), 
where coefficients C u and C,j obey the condition 2 + 2C% = 1, so that p = 

\H 2 + \H 2 + \H 2 + \H 2 = \g 0 \ 2 . 
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Figure 2 Wave packet density p at times = 0, 900, 1500, and 1800 (lattice units) for a Gaussian 
wave packet propagating through a graphene sample with randomly distributed impurities. The 
simulation is performed with impurity percentage C = 0.5% and impurity potential V = 50 MeV. 


A grid size of N z x N y = 1024 2 elements is used and the initial wave packet spread 
is set at a = 48, a fully relativistic particle (m = 0) is considered. 

In these simulations, we impose g = 0 and g = 1000 and vary the initial energy of 
the wave packet k = k z = k y in order to inspect the effect of this parameter on the 
wave packet separation, which, in turn, informs on the effective mass acquired by 
the up and down propagating modes. 

In Fig. 3, the wave function density at time t = 200 for k = 0.004, 0.04 and 0.4 is 
shown for g = 0 and g = 1000, respectively. The figure shows that sufficient energy, 
k > 0.004, is needed to observe the splitting of the wavepacket. The effects of non¬ 
linear interactions, fringes and distortions, are also well visible in the right column 
of Figure 3, A quantitative analysis in the one-dimensional case led to satisfactory 
agreement with asymptotic solutions for the dynamic mass as a function of the 
interaction strength g. A similar analysis in two spatial dimensions remains to be 
developed. 

9 Summary and outlook 

Summarizing, we have reviewed discretizations of the Dirac equation and described 
their mapping into QW’s. These relations may allow the solution of the Dirac equa¬ 
tion on quantum computers. In the first part, a general argument is given, using 
the operator splitting method. Then, the QLB scheme is studied within the same 
































16 


9 SUMMARY AND OUTLOOK 


g = o 

g = 1000 

O 

o 

O 

o 

o 

% 

o 

o 

o 

o 


Figure 3 Wave packet density p = IV’iI 2 + 1^21 2 + \^ 3 1 2 + 1^41 2 for the scheme solving Dirac 
equation with NJL interaction. These snapshots are taken at time t = 200 for k = 0.004, 0.04 and 
0.4 (from top to bottom). The figure shows how the initial energy affects the separation 
phenomenon. In particular, at low energy, k = 0.004, no splitting of the wavepackest is observed. 
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perspective and a similar relation is found. We have also shown that a similar struc¬ 
ture remains in curved space, using a scheme based on a finite volume formulation, 
with the important caveat that the exact nature of the streaming operator, typical 
of QLB, is no longer preserved. Rather, one sees the appearance of the residency 
matrix, which characterizes the fraction of spinor which is left in the cell after 
one step in the time evolution. This scheme, along with its generalization to many 
dimensions, will be studied in future work. 
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